################################### Packages ##########################################
library(tidyverse)
library(fixest)
library(stargazer)

################################### Set Working Directory ##########################################

# Getting the path of your current file
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

################################### functions ##########################################
# IHS Transformation
ihs=function(x) log(x+sqrt(1+x^2)) # to reverse ihs use sinh
# Demeans and IHS Transforms
idm=function(x) ihs(x)-mean(ihs(x),na.rm=T)
# Demeans
dm=function(x) x-mean(x,na.rm=T)



################################### Data ##########################################
### main data
up = read.csv("../Data/data_main_final.csv",stringsAsFactors = F)%>%
  ### The following codes IHS transforms the variables and demeans them at the stock level. The demeaning is only relevant for the interacted variables
  group_by(stock_id)%>% # for demeaning at the stock level
  arrange(year)%>% # for calculating the lagged variables
  mutate(
    ffmsy = ihs(ffmsy), # Use only IHS because of the simulation (below). The demeaning is only important for the interacted variables, for the non-iteracted variables nothing changes
    gdp=idm(gdp_const),
    gdp2=idm(gdp_const)^2,
    tfp=idm(rtfpna),
    fao_ag_credit = idm(fao_agriculture_credit),
    fao_to_credit = idm(fao_total_credit),
    credit=idm(domestic_credit),
    credit_instrument = idm(credit_instrument),
    interest_factor=dm(1/(1 + lending_interest/100)),
    interest_factor = ifelse(inflation_cp<=30,interest_factor,NA), # This removes observations during hyper inflation
    inflation_factor=dm(1/(1+inflation_cp/100)),
    inflation_factor=ifelse(inflation_cp<=30,inflation_factor,NA), # This removes observations during hyper inflation
    bbmsy1 = idm(lag(bbmsy,1))
    )%>%
  as.data.frame()


### ram data
ram = read.csv("../Data/data_ram_final.csv",stringsAsFactors = F)%>%
  ### The following codes IHS transforms the variables and demeans them at the stock level. The demeaning is only relevant for the interacted variables
  group_by(country, stock_id, scientific_name)%>%
  arrange(year)%>% # for calculating the lagged variables
  mutate(ffmsy = ihs(ffmsy), # Fishing mortality
         gdp=idm(gdp_const), # GDP (World Bank)
         gdp2=idm(gdp_const)^2, 
         tfp = idm(rtfpna), # Total factor productivity (Penn)
         fao_ag_credit = idm(fao_agriculture_credit), # Credit to agriculture, fisheries, forestry from the FAO
         fao_to_credit = idm(fao_total_credit), # Total credit from the FAO
         credit=idm(domestic_credit), # Domestic credit to private sector (% of GDP)
         interest=idm(lending_interest), # Lending interest IMF
         inflation=idm(inflation_cp), # Inflation from World Bank
         interest_factor=dm(1/(1 + lending_interest/100)),
         inflation_factor=dm(1/(1+inflation_cp/100)),
         bbmsy1 = idm(lag(bbmsy,1)), # Fish stocks 
         species_cat_id =as.numeric(as.factor(isscaap_group)) # Species group
         )%>%
  as.data.frame()


############## Regression for Table 2: Resource extraction and property rights ###########################
# 1. EEZ only
pr1 = feols(ffmsy ~ pr + credit + gdp + gdp2 |
              stock_id + country + year,
            cluster =  ~ country + year + scientific_name,
            data = up%>%
              mutate(pr = eez))

summary(pr1)

# 2. Main
pr2 = feols(ffmsy ~ pr + credit + gdp + gdp2 |
              stock_id + country + year,
            cluster =  ~ country + year + scientific_name,
            data = up%>%
              mutate(pr = eez*range^mobility))


summary(pr2)


# 3. Country -Year FE
pr3 = feols(ffmsy ~ pr + credit + gdp + gdp2 |
              stock_id + country^year,
            cluster =  ~ country + year + scientific_name,
            data = up%>%
              mutate(pr = eez*range^mobility))


summary(pr3)


# 4. EEZ x ROL
pr4 = feols(ffmsy ~ pr + credit + gdp + gdp2 |
              stock_id + country + year,
            cluster =  ~ country + year + scientific_name,
            data = up%>%
              mutate(pr =eez*wbgi_rle_stine))

summary(pr4)

# 5. EEZ x PRP
pr5 = feols(ffmsy ~ pr + credit + gdp + gdp2 |
              stock_id + country + year,
            cluster =  ~ country + year + scientific_name,
            data = up%>%
              mutate(pr =eez*prp_prp_stine))

summary(pr5)

# 6. EEZ x PRP
pr6 = feols(ffmsy ~ pr + credit + gdp + gdp2 |
              stock_id + country + year,
            cluster =  ~ country + year + scientific_name,
            data = up%>%
              mutate(pr = eez*wbgi_rle_stine*range^mobility))

summary(pr6)

# 7. Interaction of property rights and the rule of law index
pr7 = feols(ffmsy ~ rm +  rm:rol + rol + credit + gdp + gdp2 | 
              stock_id + country + year,
            cluster =  ~ country + year + scientific_name,
            data = up%>%
              mutate(rm = dm(eez*range^mobility), # This code demeans the variables across the whole data set (see notes in the main text)
                     rol = dm(wbgi_rle_stine)))

summary(pr7)


##### prints the results

etable(pr1,pr2,pr3,pr4,pr5,pr6,pr7,
       digits = 3,
       digits.stats = 3,
       drop = c("credit","gdp","gdp2" ),
       dict=c("pr" ="Property rights",
              "rm" = "Property rights",
              "wbgi_rle_stine"="ROL",
              "rol" = "ROL",
              "rm:wbgi_rle_stine" = "Property rights",
              "stock_id" = "Stock FE",
              "country" = "Country FE",
              "year" ="Year FE",
              "Year FE (Country FE)" = "Country trend" ),
       headers = c("EEZ","EEZ$\\times$R$^M$","EEZ$\\times$R$^M$","EEZ$\\times$ROL","EEZ$\\times$PRP","EEZ$\\times$R$^M$$\\times$ROL","EEZ$\\times$R$^M$"),
       file = "../Results/pr_final.tex",
       cluster =  ~ country + year + scientific_name,
       depvar = F,
       fixef_sizes =T,
       fixef_sizes.simplify = F,
       replace =T)



############## Regressions for Table 3: Resource extraction, property rights and credit market development #########################

### 1. Baseline without interaction term
c1 = feols(ffmsy ~ credit + pr 
           + gdp + gdp2 
           | stock_id + country + year,
           cluster = ~ country + year + scientific_name,
           data = up)

summary(c1)



### 2. Main
c2 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up)

summary(c2)



### 3. trend
c3 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year + country[year],
           cluster = ~ country + year + scientific_name, 
           data = up)

summary(c3)


#### 4. IV (the instrument starts in 1977)

c4 = feols(ffmsy ~  pr 
            + gdp + gdp2 
            | stock_id + country + year
            |credit + credit:pr ~ credit_instrument + credit_instrument:pr ,
            cluster = ~ country + year + scientific_name, 
            data = up) # Note that the credit instrument data are only available after 1977, so many observations drop out


summary(c4,stage = c(1,2))
fitstat(c4,type = "ivwald")



#### 5. OECD
# These are all OECD countries but note that the land-locked countries are not in the data (Austria, Czech Republic, Hungary,...)
oecd = c("Australia" , "Austria" , "Belgium" , "Canada" , "Chile" , "Colombia" , 
         "Costa Rica" , "Czech Republic" , "Denmark" , "Estonia" , "Finland" , "France" , "Germany" , 
         "Greece" , "Hungary" , "Iceland" , "Ireland" , "Israel" , "Italy" , "Japan"  ,
         "Latvia" , "Lithuania" , "Luxembourg" , "Mexico" , "Netherlands" , "New Zealand" , "Norway" ,
         "Poland" , "Portugal" , "Slovakia" , "Slovenia" , "South Korea" , "Spain" , "Sweden" , "Switzerland" , "Turkey" , 
         "United Kingdom" , "United States")



c5 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up%>%filter(country%in%oecd))

summary(c5)


##### 6. Balanced
years = 1960:2012
# number of years
ny = length(years)

# credit
c6 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up%>% # This code restricts the stocks to those with complete data. 
             filter(year%in%years)%>%
             group_by(stock_id,country)%>%
             filter(n()==ny) # this removes all stocks with less obsevrations than the number of years. Note that all NA values for the ffmsy variable were removed earlier
           )

summary(c6)

##### 7. Top10
### Orders species by total catch value between 2002 and 2012. We restrict this to this period to reduce the influence of random fluctuations (less years) and the influence of time series length (more years)
top = up%>%
  filter(year%in%2002:2012)%>%
  group_by(scientific_name)%>%
  summarise(value = sum(catch*price))%>% # Sums the catch values over the specified time period
  arrange(desc(value))


### selects top 10% of species
# selects all species names
spe = top$scientific_name
# number of species in the top 10 % (rounding is necessary to reduce the number to integers)
p10 = round(length(spe)*0.1,0) 
# top 10 % species
top10 = spe[1:p10]

# regression
c7 = feols(ffmsy ~ credit + pr + credit:pr 
            + gdp + gdp2 
            | stock_id + country + year,
            cluster = ~ country + year + scientific_name, 
            data = up%>%filter(scientific_name%in%top10))

summary(c7)




#### 8. Ram
c8 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, data = ram)

summary(c8)


#### 9. Ag. Credit

c9 = feols(ffmsy ~ credit + pr + credit:pr 
           + gdp + gdp2 + fao_to_credit 
           | stock_id + country + year ,
           cluster = ~ country + year + scientific_name, 
           data = up%>%mutate(credit = fao_ag_credit))

summary(c9)

###### print results
etable(c1,c2,c3,c4,c5,c6,c7,c8,c9,
       digits = 3,
       digits.stats = 3,
       drop = c("gdp","gdp2","unemployment_nat","rtfpna","fao_to_credit"),
       dict=c("pr" ="Property rights",
              "credit" = "Credit",
              "stock_id" = "Stock FE",
              "country" = "Country FE",
              "year" ="Year FE",
              "credit x pr" = "Credit $\\times$ Property rights",
              "credit:pr" = "Credit $\\times$ Property rights",
              "Year FE (Country FE)" = "Country trend" ),
       replace =T,
       order = c("!times","Property rights","Credit"),
       headers = c("Baseline","Main","Trend","IV","OECD","Balanced","Top10","RAM","Ag. Credit"),
       depvar = F,
       drop.section = c("slopes"),
       fixef_sizes =T,
       file = "../Results/main_final.tex")



######################### Regressions for Table A.1: First stage ###############################
# the main specification is c4
# here is a specification without the interaction term
# this is 
c4fs = feols(ffmsy ~  pr 
             + gdp + gdp2 
             | stock_id + country + year
             |credit  ~ credit_instrument ,
             cluster = ~ country + year + scientific_name, 
             data = up%>%
               filter(year>=1977))

summary(c4fs, stage = 1)
fitstat(c4fs,type = "ivwald")

##### Print the results
etable(c4fs,c4,
       stage = 1,
       fitstat =  c("ivwald1","n"),
       digits = 3,
       digits.stats = 3,
       drop = c("gdp","gdp2","rtfpna","fao_to_credit","pr"),
       dict=c("pr" ="Property rights",
              "credit_instrument" = "Credit instrument",
              "credit" = "Credit",
              "stock_id" = "Stock FE",
              "country" = "Country FE",
              "year" ="Year FE",
              "credit x pr" = "Credit $\\times$ Property rights",
              "credit:pr" = "Credit $\\times$ Property rights",
              "Year FE (Country FE)" = "Country trend" ),
       replace =T,
       order = c("!times","Credit"),
       headers = c("Separate","Joint","Joint"),
       depvar = T,
       fixef_sizes =T,
       file = "../Results/fs_final.tex")

###################### Bonferoni and Holm corrections for Appendix "A.11 Holm and Bonferroni p-values for the main results"
### Correction for Table A.4: Holm corrected p-values
c1h = c(round(p.adjust(c1$coeftable[c(1,2),4],method= "holm"),5),NA)
c2h = round(p.adjust(c2$coeftable[c(1,2,5),4],method= "holm"),5)
c3h = round(p.adjust(c3$coeftable[c(1,2,5),4],method= "holm"),5)
c4h = round(p.adjust(c4$coeftable[c(1,3,2),4],method= "holm"),5)
c5h = round(p.adjust(c5$coeftable[c(1,2,5),4],method= "holm"),5)
c6h = round(p.adjust(c6$coeftable[c(1,2,5),4],method= "holm"),5)
c7h = round(p.adjust(c7$coeftable[c(1,2,5),4],method= "holm"),5)
c8h = round(p.adjust(c8$coeftable[c(1,2,5),4],method= "holm"),5)
c9h = round(p.adjust(c9$coeftable[c(1,2,6),4],method= "holm"),5)

h =  tibble(c1h,c2h,c3h,c4h,c5h,c6h,c7h,c8h,c9h)%>%
  as.data.frame()

row.names(h) = c("Credit","Property rights","Credit Property rights")
names(h) = c("Baseline","Main","Trend","IV","OECD","Balanced","Top10","RAM","Ag. Credit")

stargazer(h,                 # Export txt
          summary = FALSE,
          type = "latex",
          digits =2,
          out = "../Results/holm_final.tex")

### Correction for Table A.5: Bonferroni corrected p-values

c1b =c(round(p.adjust(c1$coeftable[c(1,2),4],method= "bonferroni"),5),NA)
c2b =round(p.adjust(c2$coeftable[c(1,2,5),4],method= "bonferroni"),5)
c3b =round(p.adjust(c3$coeftable[c(1,2,5),4],method= "bonferroni"),5)
c4b =round(p.adjust(c4$coeftable[c(1,3,2),4],method= "bonferroni"),5)
c5b =round(p.adjust(c5$coeftable[c(1,2,5),4],method= "bonferroni"),5)
c6b =round(p.adjust(c6$coeftable[c(1,2,5),4],method= "bonferroni"),5)
c7b =round(p.adjust(c7$coeftable[c(1,2,5),4],method= "bonferroni"),5)
c8b =round(p.adjust(c8$coeftable[c(1,2,5),4],method= "bonferroni"),5)
c9b =round(p.adjust(c9$coeftable[c(1,2,6),4],method= "bonferroni"),5)

b = tibble(c1b,c2b,c3b,c4b,c5b,c6b,c7b,c8b,c9b)%>%
  as.data.frame()

row.names(b) = c("Credit","Property rights","Credit $\\times$ Property rights")
names(b) = c("Baseline","Main","Trend","IV","OECD","Balanced","Top10","RAM","Ag. Credit")

stargazer(b,                 # Export txt
          summary = FALSE,
          digits = 2,
          type = "latex",
          out = "../Results/bonferoni_final.tex")

########################## Section 6.3. Simulations ############################################
#### The math is for the simulation
# All variables are demeaned at the level of stock i in country c, unless noted otherwise. 
# 
# pr = [0,1] are non-transformed property rights (for stock i in country c)
# p0 := pr2012 are the property rights for stock I in year 2012
# c0 := ihs(cr2012 )– ihs (mean(cr1960,…, cr2012)) . c0 is the demeaned and ihs transformed credit for stock i in country c in the year 2012. Here, ihs (mean(cr1960,…, cr2012)) is the ihs transformed mean credit level for stock i in country c across all years for which we have data for stock i. 
# Define 
# mc = ihs (mean(cr1960,…, cr2012)) This is for stock i
# Using the estimated coefficients b,a,g, c (the fixed effect) and the variables defined above gives us the predicted ihs transformed outcome for 2012
# ihs(y0) = b p0 + a  c0 + g p0 *i0 + c  
# 
# What would happen if property rights would jump to perfect property rights and credit markets would jump to US credit market levels?
# p1 = 1 i.e. perfect property rights
# c1: = ihs(C_USA2012 )– mc C_USA2012 is the credit level in the USA in 2012. This is not demeaned. We have to subtract the stock level mean credit (mc) such that the jump in credit corresponds to the jump between the credit of country c and the US (in ihs levels).
# 
# The outcome under perfect property rights and US credit markets is then given by
# ihs(y1) = b p1 + a c1 + g p1 *c1 + c
# 
# The difference can be expressed as
# ihs(y1)-ihs(y0) = b [p1-p0] + a [c1 – c0] + g [p1 * c1 - p0 * c0]
# 
# The following must, therefore, also hold
# ihs(y1) = ihs(y0) + b [p1-p0] + a [c1 – c0] + g [p1 * c1 - p0 * c0]
# Transforming this back yields
# y1 = sinh {ihs(y0) + b [p1-p0] + a [c1 – c0] + g [p1 * c1 - p0 * c0]}  equation (*)
# 
# We can rewrite c1 as
# c1= ihs(crUSA2012 )– ihs(cr2012 ) + c0 (because the both ihs(cr2012 ) terms cancel out)
# This allows us to estimate (*) without defining the new variable mc in the dataset. Equation(8) therefore becomes
# y1 = sinh {ihs(y0) + b [p1-p0] + a [ihs(crUSA2012 )– ihs(cr2012 )] + g [p1 * c1 - p0 * c0]} 



library(ggthemes)
library(cowplot)

### Note: All variables and coefficints are as described above in the math text
# extracts the coefficients from the main specification
a=round(c2$coefficients[1],2) # credit coefficient
b=round(c2$coefficients[2],2) # pr coefficient
g=round(c2$coefficients[5],2) # credit:pr coefficient



w=up%>%
  filter(year==2012 & !is.na(pr) & !is.na(credit) & !is.na(ffmsy))%>% # uses only non-missing observations for year 2012 (baseline) 
  mutate(c0=credit)%>% # Baseline demeaned credit
  mutate(credit_us = ihs(mean(domestic_credit[which(country=="United States")])))%>% 
  mutate(c1 = credit_us - ihs(domestic_credit) + credit)%>% 
  mutate(p0 = pr)%>% # pr is non-transformed
  mutate(p1=1)%>% # maximum pr levels
  mutate(f_cr=sinh(ffmsy + (a+g*p0)*(c1-c0)))%>% # property rights remain constant such that p0=p1 and so (*) becomes: ihs(y1) = ihs(y0) + a [c1 - c0] + g [p1 * c1 - p0 * c0] = ihs(y0) + (a + g p0) * (c1 - c0)
  mutate(f_pr=sinh(ffmsy + (b+g*credit)*(p1-p0)))%>% # credit remains constant such that c0 = c1 and therefore (*) becomes ihs(y1) = ihs(y0) + (b + g c0) [p1-p0] 
  mutate(f_pr_cr=sinh(ffmsy + a*(c1-c0)+b*(p1-p0)+ # This is equation (*) 
                        g*(c1*p1-c0*p0)))%>%
  mutate(ffmsy=sinh(ffmsy))%>% # sinh(ffmsy)
  as.data.frame()


########################### Main simulation #############################
# weights by catch
wa = w%>%
  mutate(weight = catch)%>%
  filter(!is.na(weight))%>%
  summarise(
    # this just counts the number of socks that are fished beyond MSY levels
    n_cr = length(f_cr[f_cr>100])/length(f_cr)*100,
    n_pr = length(f_pr[f_pr>100])/length(f_pr)*100,
    n_pr_cr = length(f_pr_cr[f_pr_cr>100])/length(f_pr_cr)*100,
    n_fmsy = length(ffmsy[ffmsy>100])/length(ffmsy)*100,
    # this calculates the weighted mean of fishing mortality using catches as weights
    f_cr = weighted.mean(f_cr,weight,na.rm=T),
    f_pr = weighted.mean(f_pr,weight,na.rm=T),
    f_pr_cr = weighted.mean( f_pr_cr,weight,na.rm=T),
    ffmsy = weighted.mean(ffmsy,weight,na.rm=T),
  )%>%
  # The next lines reformat the results for plotting
  gather(name, value, n_cr:ffmsy)%>%
  mutate(scenario = recode(name,
                           "n_cr" = "Credit",
                           "n_pr" = "Property \n Rights",
                           "n_pr_cr" = "Credit & \n Property\n Rights",
                           "n_fmsy" = "Actual",
                           "f_cr" = "Credit",
                           "f_pr" = "Property \n Rights",
                           "f_pr_cr" = "Credit & \n Property\n Rights",
                           "ffmsy" = "Actual"),
         scenario = factor(scenario, levels = c("Actual", "Credit","Property \n Rights","Credit & \n Property\n Rights")),
         outcome = c(rep("n",4),rep("f",4)))


################################### Figure 7: Simulation results ################################
# plots mean fishing mortality
a1 = wa%>%
  filter(outcome == "f")%>%
  ggplot(aes(x=scenario, y=value)) +
  geom_bar(stat="identity", position=position_dodge(),fill="slategray4")+
  theme_bw()+
  theme(text=element_text(family="serif",size=15))+
  labs(x = "",y = "Mean extraction [%]")
a1

# Plots number of overfished stocks
a2 = wa%>%
  filter(outcome == "n")%>%
  ggplot(aes(x=scenario, y=value)) +
  geom_bar(stat="identity", position=position_dodge(),fill="slategray4")+
  theme_bw()+
  theme(text=element_text( family="serif",size=15))+
  labs(x = "",y = "Unsustainable stocks [%]")

# combines plots
cowplot::plot_grid(a1,a2, labels = c("A)", "B)") ,label_fontface="plain",nrow = 1,label_fontfamily="serif",label_size = 16)

ggsave(filename = "../Figures/simulation_agg.pdf", device = "pdf",width = 20, height = 10, units ="cm")

############################# average levels #########################################
# resource extraction
round(mean(w$ffmsy),0)
round(mean(w$f_cr),0)
round(mean(w$f_pr),0)
round(mean(w$f_pr_cr),0)

# share of overfished stocks
actual = round(length(w$ffmsy[w$ffmsy>100])/length(w$ffmsy)*100,0); actual
cr = round(length(w$f_cr[w$f_cr>100])/length(w$f_cr)*100,0); cr
pr = round(length(w$f_pr[w$f_pr>100])/length(w$f_pr)*100,0); pr
cr_pr = round(length(w$f_pr_cr[w$f_pr_cr>100])/length(w$f_pr_cr)*100,0);cr_pr

# reduction of overfished stocks
round((actual - cr)/actual*100,0)
round((actual - pr)/actual*100,0)
round((actual - cr_pr)/actual*100,0)


########################################### Figure A.10: Changes in resource extraction by species groups #############################################
library(countrycode)
w = w%>%
  mutate(region = countrycode(country,origin="country.name",destination = "region"),
         # cleans species groups
         species_group = recode(species_cat, 
                                "Mussels" = "Bivalves",
                                "Scallops, pectens" = "Bivalves",
                                "Clams, cockles, arkshells" = "Bivalves",
                                "Oysters" ="Bivalves",
                                "Miscellaneous pelagic fishes" = "Other",
                                "Miscellaneous coastal fishes" = "Other",
                                "Miscellaneous demersal fishes" = "Other",
                                "Miscellaneous aquatic invertebrates" = "Other",
                                "Miscellaneous diadromous fishes"= "Other", 
                                "Abalones, winkles, conchs" = "Other",
                                "Carps, barbels and other cyprinids"  = "Other",
                                "Miscellaneous marine crustaceans" = "Other",    
                                "Horseshoe crabs and other arachnoids" = "Other",
                                "Sea-urchins and other echinoderms" = "Sea urchins",
                                "King crabs, squat-lobsters" = "Crabs",
                                "Crabs, sea-spiders" = "Crabs",
                                "Squids, cuttlefishes, octopuses" = "Cephalopods",
                                "Herrings, sardines, anchovies" = "Herrings",
                                "Cods, hakes, haddocks" = "Cods",
                                "Lobsters, spiny-rock lobsters" = "Lobsters",
                                "Shrimps, prawns" = "Shrimps",
                                "Flounders, halibuts, soles" = "Flatfish",
                                "Tunas, bonitos, billfishes" = "Tunas",
                                "Sharks, rays, chimaeras" = "Sharks",
                                "Salmons, trouts, smelts" = "Salmon"
                                ))


# 
wb = w%>%
  group_by(species_group)%>%
  summarise(
    f_cr = weighted.mean(f_cr,catch,na.rm=T),
    f_pr = weighted.mean(f_pr,catch,na.rm=T),
    f_pr_cr = weighted.mean( f_pr_cr,catch,na.rm=T),
    ffmsy = weighted.mean(ffmsy,catch,na.rm=T),
  )%>%
  mutate(# express relative to actual
    f_cr = (f_cr - ffmsy)/ffmsy*100,
    f_pr = (f_pr - ffmsy)/ffmsy*100,
    f_pr_cr = (f_pr_cr - ffmsy)/ffmsy*100
  )%>%
  gather(name, value, f_cr:ffmsy)%>%
  mutate(scenario = recode(name,
                           "f_cr" = "Credit",
                           "f_pr" = "Property \n Rights",
                           "f_pr_cr" = "Credit & \n Property\n Rights",
                           "ffmsy" = "Actual"),
         scenario = factor(scenario, levels = c("Actual", "Credit","Property \n Rights","Credit & \n Property\n Rights")),
         outcome = ifelse(grepl("n", name),"n","f"))%>%
  as.data.frame()


b1 = wb%>%
  filter(outcome == "f" & scenario == "Credit" & species_group != "Other")%>%
  ggplot(aes(x=species_group, y=value)) +
  geom_bar(stat="identity")+
  theme_bw()+
  theme(text=element_text(family="serif",size=15),
        axis.text.x = element_text(angle = 45, hjust = 1))+
  # scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  labs(x = "",y = "Change in extraction rate [%]")+
  ylim(-45,0)
b1
  
b2 = wb%>%
  filter(outcome == "f" & scenario == "Property \n Rights"& species_group != "Other")%>%
  ggplot(aes(x=species_group, y=value)) +
  geom_bar(stat="identity")+
  theme_bw()+
  theme(text=element_text(family="serif",size=15),
        axis.text.x = element_text(angle = 45, hjust = 1))+
  # scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  labs(x = "",y = "Change in extraction rate [%]")+
  ylim(-45,0)


b3 = wb%>%
  filter(outcome == "f" & scenario == "Credit & \n Property\n Rights" & species_group != "Other")%>%
  ggplot(aes(x=species_group, y=value)) +
  geom_bar(stat="identity")+
  theme_bw()+
  theme(text=element_text(family="serif",size=15),
        axis.text.x = element_text(angle = 45, hjust = 1))+
  # scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  labs(x = "",y = "Change in extraction rate [%]")+
  ylim(-45,0)



cowplot::plot_grid(b1,b2,b3, labels = c("A)", "B)","C)") ,label_fontface="plain",nrow = 1,label_fontfamily="serif",label_size = 16)

ggsave(filename = "../Figures/simulation_species-group.pdf", device = "pdf",width = 30, height = 10, units ="cm")
ggsave(filename = "../Figures/simulation_species-group.png", device = "png",width = 30, height = 10, units ="cm")

##################################### Figure A.9: Changes in resource extraction by country #########################################
library(maps)

### the following code summarizes the results by country
v = w%>%
  group_by(country)%>%
  summarise(
    ffmsy = weighted.mean(ffmsy,catch),
    f_cr = weighted.mean(f_cr,catch),
    f_pr = weighted.mean(f_pr,catch),
    f_pr_cr = weighted.mean(f_pr_cr,catch)
  )%>%
  rename(region = country)%>%
  mutate(
    f_pr_cr = (f_pr_cr - ffmsy)/ffmsy*100,
    f_pr = (f_pr - ffmsy)/ffmsy*100,
    f_cr = (f_cr - ffmsy)/ffmsy*100,
    region = recode(region,
                         "Congo - Brazzaville" = "Republic of Congo" ,
                         "Congo - Kinshasa" = "Democratic Republic of the Congo",
                         "St. Kitts & Nevis" = "Saint Kitts",    
                         "St. Lucia"  = "Saint Lucia",
                         "St. Vincent & Grenadines" = "Saint Vincent" ,
                         "Trinidad & Tobago" = "Trinidad" ,
                         "United Kingdom" = "UK",
                         "United States" ="USA" ))


# loads word shape files and merge with simulation data
world = map_data("world")%>%
  left_join(.,v)%>%
  filter(lat>-60 & lat<85)


ggplot() +
  geom_polygon(data = world, aes(x=long, y = lat, group = group, fill = f_pr_cr), color = "black",linewidth = 0.1) + 
  coord_fixed(1.3)+
  theme_void()+
  scale_fill_gradient2(
    low = "indianred4",
    mid = "indianred1",
    midpoint = -50,
    high = "white",
    na.value = "grey95",
    limits = c(-100,0)
  )+
  theme(text=element_text(size=20,  family="serif"))+
  labs(fill = "Change in extraction\n rate [%]")

ggsave(filename = "../Figures/simulation_map.pdf", device = "pdf",width = 30, height = 15, units ="cm")
ggsave(filename = "../Figures/simulation_map.png", device = "png",width = 30, height = 15, units ="cm")


